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Abstract. In this paper we calcnlate the potential sensitivity of the CUORE detector to ax¬ 
ions prodnced in the Snn throngh the Primakoff process and detected by the inverse coherent 
Bragg-Primakoff process. The conversion rate is calcnlated rising density fnnctional theory 
for the electron density and realistic expectations for the energy resolntion and backgronnd 
of CUORE. Monte Carlo calcnlations for 5 yx741 kg=3705 kg y of exposnre are analyzed 
nsing time correlation of individnal events with the theoretical time-dependent connting rate 
and lead to an expected limit on the axion-photon conpling ga-y-y < 3.83 x 10“^*^ GeV~^ for 
axion masses less than 100 eV. 
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1 Introduction 

The CP-violating term in Quantum Chromodynamics(QCD) implies that the neutron electric- 
dipole moment should be a factor of 10^^ larger than the experimental upper bound [1], This 
CP problem was solved dynamically by Peccei and Quinn [2, 3] through the introduction of 
a global U{l)pQ symmetry. Weinberg [4] and Wilczek [5] showed that the axion, which is a 
Nambu-Goldstone boson, is produced as a consequence of the spontaneously breaking of the 
U{l)pQ symmetry. Axions, or more generally axion-like particles(ALPs), are pseudo-scalar 
bosons that can couple with the electromagnetic field or directly with leptons or quarks. The 
prospect that axions can be a candidate for dark matter in the universe [6-9] has motivated 
many experimental searches [10-21] and theoretical investigations [22, 23]. A detailed review 
of axion physics, astrophysical bounds on the mass and coupling of axions to conventional 
particles and the cosmological role of axions has been given by Raffelt [24]. 

Sikivie [25] proposed the detection of axions produced in the Sun by the Primakoff 
process shown in figure 1(a) using a magnetic helioscope. A detailed calculation of the solar 
axion flux was carried out by van Bibber et al. [26]. Buchmiiller and Hoogeveen [27] proposed 
using the Primakoff process to produce axions from an intense X-ray source by coherent Bragg 
conversion in a single crystal and the inverse coherent Bragg-Primakoff process to detect 
the axions. Paschos and Zioutas [28] proposed detecting solar axions using inverse Bragg- 
Primakoff effect, figure 1(b). GUORE(Gryogenic Underground Observatory for Rare Events) 
[29, 30] is a very low background low temperature bolometric detector that is designed to 
search for neutrinoless double beta decay(Oz//3/3). It can be sensitive enough to search for 
dark matter and solar axions. In this paper we calculate the expected sensitivity of GUORE 
to the upper bound for the coupling of axions to photons via the coherent Primakoff process 
in Te 02 single crystals. 


2 Theoretical Counting Rates 


The Lagrangian density for the electromagnetic field A coupled to the axion field </> is(in 
natural units h = c= 1) 


^=\ \ {{dtAf - 


( 2 . 1 ) 
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Figure 1. (a) An axion is produced in the solar core by the Primakoff effect: A photon couples to a 
virtual photon in the Coulomb field of the nucleus, (b) An axion couples to a charge in the detector 
via a virtual photon in the Coulomb field of the crystal producing a photon by the inverse Primakoff 
effect [12]. 


where A is a vector potential, 1/M = is the axion-photon coupling constant and cj) is 
the axion field. 

The matrix element for a conversion of an axion with momentum p and energy Ea to a 
photon with momentum k energy and polarization e can be written as: 


M 


{ke;Q\'Hint\^]p) 


1 y/a e ■ {p X k) 
{p-kf 


p{p - k)5{Ea 


E,) 


( 2 . 2 ) 


where V is the volume considered and p is the Fourier transform of the charge density dis¬ 
tribution in units of the fundamental charge e. In a periodic lattice p{p — k) vanishes unless 
p — k = G, a reciprocal lattice vector, which means that the momentum transfer q must be 
equal to G = 27r(^, j), where a, b and c are lattice constants, h, k, 1 are integers. The 

Te02 crystal used by CUORE has a tetragonal symmetry with space group 114(422) and can 
be treated as a distorted rutile structure with a = b = 4.8088 A and c = 7.6038 A [31]. 
We calculated the charge density distribution p{r) for Te02 crystal based density functional 
theory [32, 33] by utilizing the WIEN2k package [34]. p{G) is calculated by 


p{G)= f p(r)e-*^'MV 
Jv 

Table. 1 shows the comparisons between two different calculation methods, density func¬ 
tional theory(DFT) and screening length [22]. It can be seen that \p{G)\‘^ given by the screen¬ 
ing length is bigger than that given by DFT, which means that the screening length method 
gives us a higher conversion rate and thus a more strict bound on the coupling constant g'a-Y^. 
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Table 1. Selected reciprocal lattice vectors that contribute to the inverse Primakoff conversion of 
solar axions in Te 02 - 


{KKiy 

d{Kr 

^o(keV)^ 

mult^ 

\pr^\Gn 

\p^c\GY\^ 

(1,1,1) 

3.0889 

2.01 

8 

72.53 

118.98 

(2,2,1) 

1.6536 

3.75 

8 

505.88 

1021.16 

(1,2,3) 

1.6265 

3.81 

16 

220.69 

663.13 

(4,2,0) 

1.0723 

5.78 

8 

10620.90 

29125.90 

(3,2,3) 

1.1737 

5.28 

16 

481.43 

1988.35 

(3,4,1) 

0.9513 

6.52 

16 

1280.11 

3818.14 

(3,3,3) 

1.0296 

6.02 

8 

435.17 

2787.61 

(2,3,4) 

1.0842 

5.72 

16 

363.21 

114.08 

(1,5,2) 

0.9121 

6.80 

16 

24.71 

80.38 

(4,2,4) 

0.9304 

6.66 

16 

9969.21 

30530.90 

(2,3,5) 

0.9944 

7.39 

16 

549.24 

2118.81 

(5,4,1) 

0.7452 

8.32 

16 

1641.36 

5967.68 

(5,3,3) 

0.7811 

7.94 

16 

2411.45 

5240.20 

(1,3,6) 

0.9635 

6.43 

16 

15.98 

30.00 

(3,4,5) 

0.8076 

7.68 

16 

3339.10 

3984.78 

(1,5,5) 

0.7964 

7.78 

16 

2271.45 

4415.45 

(6,1,4) 

0.7265 

8.53 

16 

243.03 

117.57 

(3,3,6) 

0.8377 

7.40 

8 

47.86 

273.31 

(7,2,2) 

0.6487 

9.56 

16 

13981.86 

23033.10 

(6,0,5) 

0.7051 

8.79 

8 

984.16 

6131.84 

(2,3,7) 

0.8335 

7.44 

16 

1540.88 

2226.21 

(1,1,8) 

0.9021 

6.87 

8 

11543.29 

36814.40 

(4,4,6) 

0.7012 

8.84 

8 

238.77 

855.70 

(7,2,4) 

0.6514 

9.98 

16 

765.01 

622.35 

(6,1,6) 

0.6665 

9.30 

16 

9789.80 

27343.1 

(7,0,5) 

0.6230 

9.95 

8 

1261.85 

8147.2 

(6,3,6) 

0.6203 

9.99 

16 

14051.54 

25538.9 


Integrating over all photon final states, the lowest-order of the conversion rate becomes 


T{p) = 


AcTr'^aNc f ^ \ 


VcV \McJ 




(2.3) 


where Vc is the volume of the conventional unit cell, Nc is the number of unit cells. The cross 
section is related to the conversion rate by T{p) = ^{p)a{p) and the flux of a single axion is 

^The integers {h, k, 1) are the components of reciprocal lattice vectors G = 27r((j, t). 

is the distance between Bragg planes for a given G and d = 

^Eo is the minimum energy for which a zero rest mass particle can Bragg scatter with momentum transfer 
G. 

‘^mult is the multiplicity, or the number of reciprocal lattice vectors in each family of planes. 

is the Fourier transform of the charge density distribution p calculated by WIEN2k based on 
density functional theory. 

®|pf^(G)^| is Fourier transform of the charge density distribution p based on screening length [22]. 
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Figure 2. Solar Axion Flux with A = l(^a 77 = 10 ^GeV 


where Va is the speed of the axion. For very light axions, Va ^ c, yielding 




Pc{G) 


G2 


\p X G\ 

G2 




(2.4) 


The flux of axions from the Sun has been calculated by van Bibber et al. [26] and can be well 
approximated by the empirical form 




(2.5) 


where 


A = ( 


108Ger,4 


Mc2 


)" = (<?a77 X lO^GeVy 


( 2 . 6 ) 


is a dimensionless parameter which uses M = l^f'GeV as a benchmark, <l>o = 5.95 x 
cm“2s“i and {p{E/Eq) = exp^/^o)-i ' When helium and metal diffusion are included, 
the core temperature of the solar model will be changed a little bit. To take into account 
this small change, we use the adjusted value of Eq = 1.103 keV [22]. Figure 2 shows the 
solar axion flux due to the Primakoff process with A = 1. For convenience, a good order of 
magnitude for the counting rate can be obtained from the combination of factors 


dNo 

dE 




(2.7) 


4 











where Na is the Avogadro’s number. The number of unit cells can be expressed in terms of 
the mass of the detector m, and molar mass of the unit cell, 

m 

Nc = -Na 
l-^c 

Coherent conversion of axions to photons is possible when the energy of the axion and direction 
to the Sun, p, satishes the Bragg condition. 

Taking into account the fact that the detector has a certain energy resolution, we replace 
the delta function in eq. (2.4) with a Gaussian function ITa with the same full width at half 
maximum(FWHM) as the detector, and hnally we have the conversion rate 


dN 

~dE 


dNo dvr^a 

=mnc—— - 

dE pcVc 


G 


\p X Gp 
G6 


^[E{p,G)/Eo]Wa[E - E{p,G)] 


( 2 . 8 ) 


CUORE will have a characteristic low-energy resolution with FWHM=0.73 keV at 4.7 keV 
and a low background counting rate [35]. Finally, we integrate the total counting rate over a 
range of energies of width AF1=0.5 keV, 

fE'+AE Jjy 

R{p,E) = J^^ —{p,E')dE' (2.9) 

Figure 3 shows the calculated counting rate as a function of time over a single day for 
several energy intervals. One way to understand the time-dependent counting rate is that at 
any instant there might be one or more reciprocal lattice vectors satisfying the Bragg condi¬ 
tion. If one considers the contribution to the counting rate of a single G, one can imagine 
isodetection contours projected on the celestial sphere. Figure 4 shows the isodetection con¬ 
tours for axions with energies from 2.5 keV to 6.5 keV for G = 27r(^, ^) in steps of 0.5 keV. 

The energy bin width is chosen to be slightly bigger than the resolution of the detector. The 
cross sign at the center is the projection of G and the dotted trojectory represents a typical 
path of the Sun through that region. To give the reader some quantitative feeling for the 
angular size of the isodetection rings, the outermost ring at 6.5 keV has an angular radius 
of 72° and the ring for 6.0 keV has a radius of 70.5°, so the outermost annulus is 1.5° wide. 
The counting rate in the energy bin 6.0 — 6.5 keV will rise when the Sun passes through this 
annulus, which takes about six minutes because the Sun moves 0.25°/min. Then the next 
annulus with energy 5.5 — 6.0 keV will go up and the counting rate in the previous annulus 
will drop. 

From Figure 3 one can see that the sharpness and complexity of the counting rate 
increases with energy. This is a geometric effect because, as is clear from Figure 4, the Sun 
spends less time in the annuli with higher energy than those with lower energy and because 
there are many more reciprocal vectors available that satisfy the Bragg condition. The daily 
temporal pattern is dependent on how the Sun pass through the annuli. Figure 4 shows the 
Sun passing along a diameter and through all eight rings, so one could expect the peaks to 
be seen in all energy bins in a symmetrical pattern. On the other hand the Sun might pass 
at a grazing trojectory that crosses the outer rings without going through the inner ones, so 
the patterns shown in one energy bin may not be seen in other energy bins. 
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Figure 3. Expected counting rates Rs{E, t) of photons produced by the inverse Primakoff conversion 
of solar axions satisfied by the Bragg condition in the CUORE detector, which is located at the Lab- 
oratori Nazionali del Gran Sasso(LNGS) in central Italy(42°28'N 13°33'E). The rates were calculated 
for ^077 = ^/M = 10“® GeV~^ . 

3 Time Correlation Method 

The distinct variation of the connting rates as a fnnction of time snggests that the detection 
of solar axions can be analyzed nsing the time correlation method^: 

N 

X = '^W{ti) X n{ti) (3.1) 

i=l 

where n{ti) is a Poisson-distribnted random variable at time tj, W(ti) is the weighting fnnction 
at time ti, which is snbject to the following two constraints 

W{t)dt = 0 (3.2a) 

W‘^{t)dt < oo, (3.2b) 


^Here we suppress the index referencing to a particular energy bin for simplicity 
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Figure 4. Isodetection contours projected on the celestial sphere for the G = 27r(-, - , plane. 
The cross sign at the center is the projection of the normal to the (1,1,1) plane. The dotted trojectory 
represents the path of the Sun through that region. 


The expected number of counts in a time interval Ati is: 

{n{ti)) = {Rbg + ^Rsi'ti))Ati 


(3.3) 


where Rbg is the background counting rate, Rsit) is the theoretically expected counting 
rate of solar axions at ga^y'y = 10“® GeV~^ and A is defined in eq. (2.6). Note that n(tj) is 
essentially equal to 1 or 0 if Ati is very small. Then the average of y is simplified to 


(X(A)) = A / W{t)Rs{t)dt 
Jo 


(3.4) 


The variance of x becomes 


(Ax(A)) 2 = (x^) - (x)^ 

rT 


= Rbg W\t)dt + X 
Jo 


W^{t)Rsit)dt 


(3.5) 


Note that the rmcorrelated events are canceled out and each event follows Poisson statistics: 

{n^{ti)) - = {n{ti)) (3.6) 

The second term in eq. (3.5) is negligible compared with the first term when A is small, so 


Ax' = Rbg / W\t)dt 


(3.7) 
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The number of events in each time interval is statistically independent, so the probability 
distribution of y given a weighting function W(t) is 


P{x\w) = ^ 


(3.8) 


The delta function can be represented by its Fourier transform, 

N 


P{x\W)= r — 

J-^ 27r \ / 


(3.9) 


Expanding the average about w = 0 and keeping the first two terms gives, by eq. (3.4) and 
eq. (3.5) 


/■oo ^ 

p(xiH-)=y_^ 


exp 


-*w(x- (x)) - 


duj 


\/27rAx^ 


exp 


(x - (x) ) 
2Ax2 


21 


(3.10) 


which shows that the probability distribution for x given the weighting function W{t) is a 
Gaussian. This is an example of the Central Limit Theorem, which states that data which are 
affected by many small and unrelated random effects are approximately normally distributed. 
We want to choose W{t) to maximize (x) subject to the constraints in eq. (3.2). Using the 
method of Lagrange multipliers we want to maximize 


^ = {x) — d'l [ W{t)dt — ^2^ [ W‘^{t)dt 

Jo Jo 

with respect to W (t), which gives 

Rs{t) - /ii - 2 ^ 2 AlT(t) = 0 
where and fj ,2 are multipliers. So 

can be determined by using the constraint that W(t}dt = 0 


(3.11) 


(3.12) 


(3.13) 


/ (Rsit) - fii)dt = 0^ = Rg 

Jo 

where Rg is the average of Rs{t) over the time considered. The second Lagrange multiplier 
determines the norm of the weighting function. It is convenient to choose /U 2 = ^ so that 
the “best” weighting function is 


W{t) = Rgit) - Rg{t) 


The log likelihood function for A is 


L(A) oc — 


{x-XfW{t){R8-Rg)dty 

2Rbg W^{t)dt 


(3.14) 


(3.15) 
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The most probable value for A is 


ix) 


(3.16) 


A = 


The width of the likelihood function is 

RBG!W\t)dt 

{jW{Rs-Rs)dtf 

Taking the variational derivative with respect to W then gives 

9AA2 W f WjRs - R8)dt - f W^(t)dt(R8 - Rg) 
(fW(R8 - R8)dt)^ 

which vanishes if IT = R8 — Rs- With this choice of weighting function 


(3.17) 


(3.18) 


\/ fW^t)dt 

if we choose the weighting function to be eq. (3.14), not only is (x) maximized but also 
AA is minimized. The generalization to the case with several independent energy bins is 
straightforward; eq. (3.16) becomes 


Efc(Afc) 

Ekfwiit)dt 


and eq. (3.19) becomes 

AA = 

where k is the index for the energy bins. 


Rbg 




4 Monte Carlo Simulation 

In a real experiment the total counting rate is given by 

R = Rbg + 


(3.20) 


(3.21) 


(4.1) 


The background counting rate will dominate when A is small. We can evaluate the sensitivity 
of the time correlation method by generating pseudo-random data with a given value of A 
and then evaluating A and AA from eq. (3.15), eq. (3.20) and eq. (3.21). 


4.1 Generation of Pseudo-data 

Let Po(Lto) be the probability that no event occurs from to to t. For a time-dependent 
counting rate i?(t), 

Po(t + At, to) = Po(t, to) X (1 - R{t) X At) 

which leads to 

Po(i,to) = (4.2) 
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The probability that the first event takes place at time ti in a small time interval At is 

P^{ti) = e~ R{ti)At 


In order to generate a sequence of events with the proper probability distribution, a 
random number M uniformly distributed on the interval 0 < .^ < 1 is used to determine the 
time t by solving 


F{t) = I - e-fo 


(4.3) 


Note that F(t) also lies in the interval 0 < F{t) < 1. The probability of choosing a random 
number ^ in a small interval AS^ is equal to the probability of hnding an event at ti in a 
small time interval At 

A ^ dF 

AM = -^At 

dt 


Differentiating F{t) gives 


AM = Po{ti,0)R{ti)At = Pi{ti)At 


Hence the time of the hrst event will be distributed correctly by solving F{ti) = M. At this 
point the time is reset and the procedure repeated until we reach the end of the simulation. 
The times {ti, i = 1, - ■ ■ ,N} are used to calculate x = ^(^*)) from this we extract 

A and AA as described above. 


5 Conclusions 

We have carried out Monte Carlo calculation using the mass, energy resolution and realistic 
background for the CHORE detector operating for 5 years. The results show that the CHORE 
detector with 741 kg Te02 in operation for 5 years can set an upper bound on A of 

A < 2.15 X 10“® (5.1) 

which is equivalent to an upper limit on the axion-photon coupling constant < 3.83 x 
10“^*^ at 95% conhdence level. To illustrate the resolving power of the time correlation 
method, in hve years with = 3.83 x 10“^® there are approximately 600 events due to 
axion conversion and 5.5 x 10^ background events. 

Figure 5 is an exclusion plot comparing this calculation with the best limits set by CAST 
[17, 18, 21] on the gayy-ma plane. The lightly shaded area and dotdashed line correspond to 
various theoretical axion models [36-39]. Our predicted bound is comparable to the newest 
CAST results for axions with mass less than 1.2 eV [21] and will improve the bound for axion 
masses in the range 1 eV < rua < 100 eV indicated by the darker shaded region(green in 
color). 

Recently, the International Axion Observatory (lAXO), a new generation axion helio¬ 
scope searching for solar axions by Primakoff conversion in a strong magnetic held, has been 
proposed (see recent work by J. K. Vogel et al. [40]). The predicted sensitivity of lAXO 
to the coupling constant Qayy is predicted to be on the order of 4 x 10“^^ GeV~^ for axion 
masses less than 0.1 eV. This is a great improvement over all current experiments, narrowing 
down search region for axions and dark matter signihcantly. However, this excluded region 
of parameter space will not reach beyond 0.2 eV. 

^The upper limit of 100 eV is somewhat arbitrary and conservative. The Bragg conversion probability is 
not very sensitive to axion masses less than 100 eV, and solar axion flux also varies very little. For axion masses 
of several hundred eV the solar axion spectrum is distorted and decoherence begins to affect the conversion 
probability. 
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Figure 5. Exclusion limits on the Qa-y-y-'nia plane. The shaded area is favored by the KSVZ [36, 37] 
and the DFSZ [38, 39[ axion models. The dotted line shows that with 3705 kg y of data, CUORE 
could exclude axions with ^£,77 > 3.83 x 10“^° GeV~^ and masses less than 100 eV. 
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